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Fullerene like cages and naonotubes of carbon and other inorganic materials are currently under 
intense study due to their possible technological applications. First principle simulations of these 
materials are computationally challenging due to large number of atoms. We have recently developed 
a fast, variational and fully analytic density functional theory (ADFT) based model that allows study 
of systems larger than those that can be studied using existing density functional models. Using 
polarized Gaussian basis sets (6-311G**) and ADFT, we optimize geometries of large fullerenes, 
fullerene-like cages and nanotubes of carbon, boron nitride, and aluminum nitride containing more 
than two thousand atoms. The calculation of C216O using nearly 39000 orbital basis functions is the 
largest calculation on any isolated molecule reported to-date at this level of theory, and it includes 
full geometry optimization. The electronic structure related properties of the inorganic cages and 
other carbon fullerenes have been studied. 
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Computer simulations are playing increasingly impor- 
tant role in our understanding about materials. Gen- 
erally, the choice of computational models that are em- 
ployed in studying the properties of materials depend 
on the property of interest and the length scale or the 
size of the system|l|. The latter is the most impor- 
tant factor in the selection of appropriate level of theory. 
Our interest is in the electronic and structural proper- 
ties of large carbon fullerenes and fullerene like cages of 
aluminum and boron nitride containing a few hundred 
atoms. At these sizes, the current toolbox of methods 
that are available include semiempirical quantum me- 
chanical models such as ZINDoHl, PM3y| methods or 
tight binding approaches 4| . More accurate description of 
electronic properties require use of more involved meth- 
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ods such as density functional (DF) based models 




The traditional quantum chemical beyond Hartree-Fock 
methods or quantum Monte Carlo are, in general, more 
accurate than the DF models. However, they are suitable 
for systems with a few tens of atoms. At present, the 
applicability of DF models is restricted to two to three 
hundred atoms depending on the schemes used to ap- 
proximate kinetic and exchange energy functionals, the 
basis sets used to expand the Kohn-Sham orbitals, the 
treatment of core electrons (use and quality of pseudopo- 
tentials), and the type of atoms in the system. We are 
working towards development of fully analytic implemen- 
tation of density functional theory ( ADFT) 0, Q. The 
computationally efficient ADFT and efficient use of the 
available point group symmetry of molecules allow us to 
optimize large inorganic and carbon fullerenes containing 
more than two thousand atoms [illlo||. 

The ADFT uses analytic atom-centered, localized 
Gaussian basis sets. These basis sets are used to to ex- 
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pand the molecular (Kohn-Sham) orbitals and the one 
body effective (Kohn-Sham) potential using variational 
and robust fitting methodology^! Q. The exchange- 
correlation part of Kohn-Sham potential is obtained us- 
ing the functional form that is based on Slater's ex- 
change functional^. For this reason, the analytic 
implementation is also called as Slater-Roothaan (SR) 
method 7] . The SR method allows an arbitrary scal- 
ing of the exchange potential around each type of atoms 
in the heteroatomic systems. These scaling factors can 
be used to parametrize the SR method. Using a suit- 
able choice of these scaling parameters, accurate total 
and atomization energies that are comparable to some 
of the most so phis ticated density functional models can 
be obtained^ LLJ, Ljj- In the following section we de- 
scribe the analytic implementation of the density func- 
tional model and the details of the SR method. 



A. Analytic formulation of the 
Gaspar-Kohn-Sham-Slater density functional model 

In the Hohcnbcrg-Kohn-Sham formulation of the den- 
s , y functional tta^fl the total electee ene^ of 
system containing N electrons and M nuclei is given by 

N 

E HKS [p] =J2< > +E ee +E xc [ PhPl ] (1) 

i 

where, the first term contains the kinetic energy operator 
and the nuclear attractive potential due to the M nuclei, 
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The second term in Eq. (JJJ represents the classical 
Coulomb interaction energy of electrons while the last 
term is the exchange-correlation energy that represents 
contributions that are quantal in origin. The Eq. Q 
is an exact expression for the total energy but practi- 
cal application require approximation to the E xc . Over 
the years numerous parameterization of different accu- 
racy and complexity have been devised and are avail- 
able in literature to model E xc . Most of them however 



have complex functional form making use of numerical 
grids necessary in implementations of the DF models. 
Number of groups have developed numerical integration 
methods for computation of integrals over the exchange- 
correlation contributions 16]. Today practically all im- 
plementations of the DF models use numerical grids to 
compute the contribution to the total energy and matrix 
elements from the exchange-correlation terms. This is 
true even if analytic basis sets such as Gaussian are used 
to express KS orbitals. However, it turns out that if one 
models the E xc according to Gaspar-Kohn-Sham-Slater 
then the contribution to total energy from this term can 
also be calculated analytically using the Gaussian basis 
sets and variational methodology [la IttI ll^] . 

The Gaspar-Kohn-Sham-Slater (GKS) exchange en- 
ergy functional is given by 
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p|(r) + pf(f) . (3) 



where, a — 2/3 is the Gaspar-Kohn-Sham value and 
a = 1 is the Slater's value. In order to calculate E xc 
analytically the one-third and two-third powers of the 
electron density are expanded in Gaussian basis sets: 



P 3 (r) 
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Here, {Ei} and {Fi} are independent Gaussian basis 
functions, while e< and /j are expansion coefficients. The 
exchange energy is then given by^J 0, list 
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where C a = -9a ^ . Thus using four LCGO basis 
sets (one for orbital expansion and three for the fitting 
the Kohn-Sham potential) the total energy is calculated 
analytically. 

Similarly, to compute the Coulomb energy, 



Eee = (P\\P) = 



p(r)p{r') 3 3 , 

— —a rd r , (7) 



we use the first robust and variational fitting methodol- 
ogy and express the charge density as a fit to a set of 



Gaussian functions, 



(8) 



Here, p(r) is the fitted density, d^ is the expansion coef- 
ficient of the charge density Gaussian basis- function G^. 
The elimination of the first order error in total energy due 
to the fit leads to the unique robust expression for the 
self-Coulomb energy [nj. The LCAO orbital coefficients 
and the vectors d, e, and f are found by constrained vari- 
ation. It is easy to obtain contribution from the first term 
in Eq. (JJJ in analytic fashion. Thus, in ADFT four sets 
of Gaussian basis are required: one for KS orbitals, three 
for the KS effective potential. This methodology was 
successfully implemented by Werepentski and Cook who 
demonstrated that noise- free forces and smooth potential 
energy can be obtained using a fully analytic (grid-free) 
implementation jl^. Q] . 



the SR method has following form 
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Here, C x = C a /a; the partitioned 3/4 power of the ex- 
change energy density, 
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where Dfj (r) is the diagonal part of the spin density ma- 
trix multiplied by the partitioning function, 



a(i) = a 



3/8 



(11) 



which contains on , the a in the Xa model for the atom on 
which the atomic orbital i is centered. The fits to powers 
of <? CT are obtained variationally from Eq. QjJj . 



COMPUTATIONAL DETAILS 



B. Slater-Roothaan method 

While the above analytic implementation is computa- 
tionally efficient its performance for the atomization of 
molecules is limited due to the limitation of the func- 
tional form adopted. We have tested its performance by 
computing atomization energies of a set of 56 molecules 
from the G2 dataset. For a = 2/3, the mean absolute er- 
ror in atomization of 56 molecules is about 16 kcal/mol. 
This can be improved to 12 kcal/mol by allowing value of 
a to change 19]. Thus the a in Eq. © can be viewed as 
a scaling parameter that scales exchange potential. The 
above model then can be modified so that each type of 
atom in the heteroatomic system has its own value of 
scaling parameter. This led to the development of the 
Slater-Roothaan (SR) method 7]. Apart from the advan- 
tage that the calculations can be performed in complete 
analytic fashion, it also allows molecules to dissociate 
c„„« tly ,„ atoned taftQ 

. The exchange energy in 



As noted earlier the analytic SR method requires four 
Gaussian basis sets. One for the orbital expansion and 
others to fit different powers of electron density, which 
we obtain from literature. We choose Pople's triple-^ 
(TZ) 6-311G** basisOJQ and the DGaussQ valence 
double-C (DZ) basis set [23 called DZVP2 for orbitals ba- 
sis sets. The s-type fitting bases are obtained by scaling 
and uncontracting the s part of the orbital basis. The 
scaling factors are 2 for the density, | for and | 

for p3 . These scaled bases are used for all s-type fitting 
bases. Ahlrichs' group has generated a RI-J basis for fit- 
ting the charge density of a valence triple- (T orbital basis 
set used in the Turbomole program |26J. The non-s 
parts of Ahlrich's fitting bases are used in combination 
with the 6-311G** orbital basis sets. We use this com- 
bination of basis sets (6-311G**/RIJ ) for boron nitride 
cages and carbon fullerenes. In combination with DZVP2 
orbital basis, we use the pd part of the A2 charge den- 
sity fitting basis. The combination DZVP2/A2 is used 
for studying aluminum nitride cages. The geometries 




FIG. 1: The optimized BN cages: (a) Two views of the 
octahedral B24N24 cage, (b) Two views of the S4 B24N24 
cage, (c) Ss B24N24 cage, (d) B28N28 cage of T symmetry, 
(e) B36N36 cage of T<j symmetry, (f) octahedral BggNg6 
cage, and (g) hemispherical cap of (8,8) BN nanotube based 

on half of Bg6Ng6 . 

of molecules were optimized using the Broyden-Fletcher- 
Goldfarb-Shanno (BFGS) algorithm Q- The forces on 
atoms are rapidly computed non-recursively using the 4-j 
generalized Gaunt coefficients 0, 0] . The atomic ener- 
gies are obtained in the highest symmetry for which the 
self-consistent solutions have integral occupation num- 
bers. The atomization energy is computed from the to- 
tal energy difference of optimized molecule and its con- 
stituent atoms. 

II. BORON AND ALUMINUM NITRIDE CAGES 

The discovery of carbon fullerene Cgo , followed by dis- 
covery of higher fullcrcncs and carbon nanotubes has led 
to intense search for hollow cage-like and tube-like struc- 
tures of other materials. In this search, boron nitride 
(BN) is probably the second most studied material after 
carbon. A number of groups have reported observation 
of BN nanotubes as well as cage like structures 

Particular relevant to this article are the series of ex- 




FIG. 2: The optimized structures of capped BN nanotubes. 
(a) and (b) are two different views of the B24N24 cage: (a) 
along the C3 axis, (b) along C4 axis, (c) B28N28 {Cak) cage 
obtained by adding a ring of 8 alternate B and N atoms, (d) 
tubular B32N32 (Ss) cage by inserting two rings of 8 alternate 
B and N atoms (see text for more details). 



TABLE I: The calculated values of binding energy per A1N 
pair, the energy gap between the highest occupied molecu- 
lar orbital and the lowest unoccupied molecular orbital, the 
vertical ionization potential (VIP), the electron affinity, and 
the energy gap obtained from the ASCF calculation for the 
optimized A1N cages. Last row gives range of values for the 
same set of BN clusters. All energies are in eV. 

Symmetry BE GAP VIP VEA A SCF 



A124N24 





10.24 


2.97 


7.05 


1.46 


5.59 


A124N24 


S 4 


10.34 


2.47 


6.84 


1.72 


5.12 


A124N24 


Ss 


10.34 


2.63 


6.79 


1.58 


5.21 


A1 28 N28 


Cih 


10.42 


2.74 


6.81 


1.59 


5.22 


A1 28 N28 


T 


10.45 


2.67 


6.84 


1.69 


5.15 


A132N32 


S s 


10.49 


2.79 


6.77 


1.61 


5.16 


AlseNse 


T rf 


10.54 


2.70 


6.73 


1.76 


4.95 


A1 48 N 4 8 


Sd 


11.09 


2.81 


6.56 


1.76 


4.8 


A1 96 N 96 





10.48 


2.18 


6.15 


2.34 


3.8 


BN range 




15 


4-5 


7-9 





7-9 
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periments by Oku and coworkers in which they detected 
BN clusters in mass spectrum [29! I^ j. These authors 
have proposed a number of cage like structures for the BN 
clusters detected in mass spectrum. Here, we report the 
electronic structure of these cages and their aluminum 
nitride (A1N) analogues. We note that while of the BN 
cages has been reported, no cage like structure of A1N 
have not yet been observed although observations of the 
A1N nanotubcs have been recently noted Q EH Q 
The optimized cage structures of BN are given in Fig. 
n Also, given are the symmetries of these cage struc- 
tures. All these cage structures have been found to be 
energetically stable with binding energy (BE) of about 
14-16 eV per pair of BN. Notable amongst these is the 
octahedral B24N24 cage that was proposed by Oku and 
coworkers as a candidate structure for one of the most 
abundant cluster in the mass spectrum. This cage is per- 
fectly round and like in Cgo fullerene where each carbon 
atom is equivalent to all other carbon atoms, a pair of 
BN in this cluster is equivalent to other pairs in the clus- 
ter. It is to be noted that the exact analogue of carbon 
fullerene is not possible for alternate boron nitride cages. 
The presence of pentagonal rings in carbon fullerenes do 
not permit full alternation of B and N atoms. Thus 
even membered rings are necessary to make alternate 
fullerenes close. The octahedral round B24N24 cage con- 
tains six square and six octagons. This structure can be 
used to form caps for (4,4) BN nanotubes 37]. However, 
unlike Ceo fullerenes, the round B24N24 octahedral cage 
is not energetically special. The other alternate B24N24 
cages with symmetry S4 and S8 are energetically nearly 
degenerate with octahedral cageQ. So it is not clear 
that which structure is likely to be observed in the ex- 
periment. 

The C4h B 2 8N 2 8 cage [Fig.^c)] can be generated from 
the base B24N24 cage by cutting the latter into two halves 
after orienting it along the C4 axis, then inserting a ring 
of eight alternate B and N atoms perpendicular to the 
axis, i.e. horizontally, and then rotating the top half by 
an eighth of a revolution. The resultant cage contains 8 
inequivalent atoms and has C^h symmetry If two rings 
of four alternate BN pair are inserted instead of one then 





FIG. 3: Optimized octahedral Alg6N96 cages: (a) Two shell 
onion-like octahedral cage with AI24N24 cage at its interior 
AI96N96 -I, (b) Fullerene-like cage AI96N96 -II. 



the resultant B32N32 cage is a tubular structure with Ss 
symmetry. The binding energy systematically increases 
by going from B24N24 to the B 2 8N 2 8 cage (0.26 eV per 
BN pair) and from B 28 N28 to B 32 N 3 2 cage (0.06 eV/BN 
pair). The successive additions of alternate BN rings 
energetically stabilizes the BN tubular cages and results 
in (4,4) BN nanotube with round caps that are based 
on octahedral B24N24 cage. Note that same (4,4) tube 
can also be generated by starting with S8 B24N24 cage 
structure. The hemispherical caps of (4,4) tube based on 
octahedral B24N24 that we have proposed have also been 
observed in a molecular dynamics study of the growth 
mechanism of BN nanotubes |39|. 

The B24N24 can be enlarged by adding hexagons. This 
leads to another round cage B96N96 of octahedral sym- 
metry. The optimized BggNgg cage is shown in Fig. ^ 
(f). Energetically the BggNg6 cage is more stable than 
the B24N24 cage. It is clearly different, from B24N24 , in 
that while being mostly round, its twelve squares stick 
out significantly, like the detonators of a sea mine. Its 
halves can form a round cap for the (8,8) BN nanotube 
(See Fig. IU(g)). 

We have reoptimized the BN cages by replacing boron 
by aluminum. We would like to point out that unlike BN 
cages which are experimentally observed, the A1N cages 
studied here are predictions. The optimized cage struc- 
tures in A1N are similar to those of BN except that they 
are larger in size due to the larger A1N bond distance than 
that of BN. The exception to this trend is the AlggNgg 
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TABLE II: The median nearest-neighbor bond distance, 
mean radius, radial standard deviation, all in Angstroms, 
for the fullerenes of this work computed using = 0.684667. 
The right-most column gives the atomization energy per atom 
(AE), in electron volts, that we compute using = 0.64190. 
Fullerene Median bond distance Mean radius AE 



C60 


1.4244 


3.5481 


-7.140 


C240 


1.4306 


7.0728 


-7.373 


C54O 


1.4264 


10.5528 


-7.431 


C96O 


1.4249 


14.0342 


-7.459 


C1500 


1.4244 


17.5225 


-7.474 


C216O 


1.4241 


21.0137 


-7.484 



Fully optimized fullerenes by analytic DFT (all electron calc. 63110" basis) 




cluster (See Fig. 0). The optimization of AI96N96 struc- 
ture starting from B96N96 cluster results in formation of 
double shell onion like structure. This onion cluster has 
AI24N24 cage at its core 40J. On the other hand, if one 
scales the B36N36 to account for larger A1N bond length 
and then replace B by Al and optimize then one does get 
fullerene like hollow cage with squares sticking out. This 
cage will be refereed to as AI96N96 -II. We find that all 
A1N cages are energetically stable with binding energy of 
about 10-11 eV per pair of A1N. However, the binding 
energy of A1N cages is less that BN cages, which have 
binding energy of 14-16 eV per BN pair. Similarly, the 
ionization potential (IP) and the energy gap between the 
energy eigenvalues of highest occupied molecular orbital 
(HOMO) and the lowest unoccupied molecular orbital is 
smaller in A1N cages. The vertical IP of BN cages is in 
the range 7-9 eV while it is about 6-7 eV in A1N cages. 
Due to the large HOMO-LUMO gap in the BN cages, 
the BN cages do not like to bind an extra electron. Con- 
sequently, the electron affinity of boron nitride cages is 
practically zero within our model. The A1N cages have 
electron affinity of about 1-2 eV. We have summarized 
the electronic structure data of A1N cages in Table I. In 
the last row of the same table contains the range of values 
for the BN cages. The HOMO-LUMO gap calculated by 
the so called A SCF method in which the first ionization 
potential is subtracted from the first electron affinity is 
given in the last row of the table. 



FIG. 4: Fully optimized structures carbon fullerenes (Basis 
set: 6-311G**/Ahlrichs) (see text for more details). 

III. CARBON FULLERENES 

Carbon fullerenes structures larger than 100 atoms 
have been studied by several groups^! I42I l^ . Most 
of these studies have used semiempirical models or tight 
biding methods or the Hartee-Fock theory plus minimal 
basis sets. Except for very recent calculation [Zl) on C240, 
fullerenes have not been studied using reasonable quality 
basis sets 0| . This is principally because of large com- 
putational cost. We have used computationally efficient 
ADFT described above to optimize geometries of several 
carbon fullerenes from Cqq to C2160 using large polarized 
Gaussian basis sets of triple zeta quality f6-311G**')[lo|. 
The ADFT code developed in our group exploits the 
icosahedral symmetry of these fullerenes in an efficient 
manner. Therefore, very large calculation on C2160 with 
about 39000 orbital basis functions, is still doable with 
modest computation resources. In order to get accurate 
geometries of larger fullerenes, we parametrize the ADFT 
to get the exact geometry of Cgo- This is accomplished 
by minimizing the mean square deviation between the 
experimental and predicted bond lengths of Cgo- This is 
possible without much difficulty as optimization of Ceo 
using ADFT takes less than 5 minutes on single proces- 
sor Linux box (Intel(R) XEON(TM) CPU 2.20GHz with 
2Gigabytes of Random Access Memory). The exact ge- 
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ometry of Ceo can be obtained for a = 0.684667. We 
use this value of a for optimizing larger carbon fullerenes 
and hope that this will give accurate estimates of their 
geometries. The C960 fullcrene can be optimized on sin- 
gle processor in 5 days. The C2160 optimization was 
performed on Linux cluster using 48 processors and took 
about 5 days. The median bond distance of optimized 
carbon fullerenes is given in Table II and optimized struc- 
tures are shown in Fig. 01 We also made an attempt to 
get accurate atomization energies using the optimized ge- 
ometries of fullerenes. For this purpose we reparametrize 
ADFT to get the exact binding energy of Ceo fullerene 
and use the a value thus determined to compute the at- 
omization energies of larger fullerenes. These values are 
also given in Table II. However, such a procedure fails 
in that the atomization energy of C240 is already lower 
than that of graphite. Thus, to get accurate estimate of 
atomization energy we need to go beyond the functional 
form that we have chosen in parameterizing the ADFT. 



Work is progress in our laboratory in this direction. 



IV. CONCLUSION 

We have presented fully analytic implementation of 
density functional theory (ADFT). It uses analytic Gaus- 
sian basis sets to varitationally express the Kohn-Sham 
molecular orbitals, electron density and the one body ef- 
fective potential of the density functional theory. The re- 
sultant formulation is computationally very efficient and 
allow for calulations on relatively large systems. It per- 
mits use of atomic number dependent potential by means 
of Slater's exchange parameters. Using the ADFT code, 
which efficiently uses the full point group symmetry of 
the molecule, we have optimized large inorganic fullerene- 
like cages and carbon fullerenes containing more than two 
thousand atoms. 
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